function [vk, vep, ck, cX, cu] = fn_solve_tba(T,H,mumu,Nc,Num_point)


h0 = H/T;
Num = Num_point;

vk = linspace(0, 3*sqrt(mumu + T*log(2*cosh(h0/2))), Num);
% vep_ori = vk.^2 - mumu - T*log(2*cosh(h0/2));
vep_ori = vk.^2 - mumu - T*log(2*cosh(H/2/T));

% F = @(x, T) -T*log(1+exp(-x/T));
% F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
%     ./(1 + fn_inveta(m, h0)));
ck = cell(1, Nc);
cX = cell(1, Nc);
cu = cell(1, Nc);

for n = 1:Nc
    ck{n} = linspace(0, 4, Num);
    cX{n} = zeros(1, Num);
    cu{n} = zeros(1, Num);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vep = vep_ori;

while 1
    
    vep0 = vep;
    ck0 = ck;
    cX0 = cX;
    
    vep = vep_ori + fn_an(2, vk, fn_FF(vep, T), vk);
    
    for n = 1:Nc
        vep = vep + T*fn_an(n, ck{n}, cX{n}, vk);
    end
    
    vnorm = zeros(1, Nc);
    for n = 1:Nc
        
        vFT = -1/T*fn_FF(vep, T);
        [vp, vr, vu] = fn_iteration(n,h0,vk, vFT, ck, cX, Num);
        ck{n} = vp;
        cX{n} = vr;
        cu{n} = vu;
        
        vnorm(n) = fn_norm_2(ck{n},cX{n},ck0{n},cX0{n});
        
    end
    
    diss1 = fn_norm(vep - vep0);
    diss2 = fn_norm(vnorm);
    fprintf('diss = %1.2e   ', diss1);
    fprintf('diss = %1.2e   \n', diss2); 
    
     if diss1<1e-8 && diss2<1e-8
        break;
    end

end

end